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Abstract 

The flow of power law fluids, which include shear thinning and shear thick- 
ening as well as Newtonian as a special case, in networks of interconnected 
elastic tubes is investigated using a residual based pore scale network model- 
ing method with the employment of newly derived formulae. Two relations 
describing the mechanical interaction between the local pressure and local 
cross sectional area in distensible tubes of elastic nature are considered in 
the derivation of these formulae. The model can be used to describe shear 
dependent flows of mainly viscous nature. The behavior of the proposed 
model is vindicated by several tests in a number of special and limiting cases 
where the results can be verified quantitatively or qualitatively. The model, 
which is the first of its kind, incorporates more than one major non-linearity 
corresponding to the fluid rheology and conduit mechanical properties, that 
is non-Newtonian effects and tube distensibility. The formulation, implemen- 
tation and performance indicate that the model enjoys certain advantages 
over the existing models such as being exact within the restricting assump- 
tions on which the model is based, easy implementation, low computational 
costs, reliability and smooth convergence. The proposed model can there- 
fore be used as an alternative to the existing Newtonian distensible models; 
moreover it stretches the capabilities of the existing modeling approaches to 
reach non-Newtonian rheologies. 

Keywords: fluid mechanics; power law fluid; elastic networks; elastic porous 
media. 


1 Introduction 

The flow of non-Newtonian fluids in networks of interconnected distensible con- 
duits occurs in many fluid dynamic systems especially the biological ones in living 
organisms. The flow of blood; which is a complex suspension with non-Newtonian 
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attributes that include deformation rate dependency, yield stress, thixotropy and 
viscoelasticity; is an obvious example [1-4], 

Although there are some previous investigations related to the non-Newtonian 
flow effects in small distensible branching flow ensembles [5-7], we are not aware 
of an attempt to deal with this issue in a direct, general and large scale approach. 
The obvious reason is the theoretical and computational difficulties associated with 
modeling such a complex phenomenon which contains multiple non-linearities re- 
lated to the rheology of the conduit wall as well as the rheology of the fluid itself plus 
fluid-structure interaction. The existing and commonly used Navier-Stokes models, 
whether the one-dimensional or the multi-dimensional ones, and their derivatives 
fall short of incorporating such effects due to their restriction to the Newtonian 
flow. Incorporating non-Newtonian effects in such models does not only require 
major model development but it also requires considerable approximations that 
normally compromise the results. 

In this paper, we investigate the flow of shear-thinning and shear-thickening 
time-independent non-Newtonian fluids of power law type in elastic networks using 
newly derived formulae [8] in conjunction with the residual based pore scale network 
modeling [9]. The Newtonian case is obtained as a special case by setting the flow 
behavior index of the power law fluid to unity. This inclusion is confirmed by the 
results obtained from specially derived formulae for the Newtonian case. 

Based on this approach, the pressure and volumetric flow rate solutions in such 
networks are obtained by implementing the characteristic flow equations for the 
particular fluid and the defining mechanical properties of the network tubes in a 
residual based pore scale network modeling code as described in detail in [9] and 
outlined in section 2. In our formulation and implementation, we use two elastic 
models to characterize the mechanical response of the tube wall that correlates the 
magnitude of the local pressure to the magnitude of the local cross sectional area. 

In modeling the network flow, we adopt the same assumptions on which the 
derivation of the analytical tube equations [8] used to characterize the flow is based, 
that is the flow in each tube is presumed incompressible, laminar, time-independent, 
fully developed at comparatively low Reynolds numbers with minimal entry and 
exit edge effects. Hence, the flow is regarded as essentially viscous with minimal 
inertial contributions. 

Furthermore, the walls of the tubes are assumed to be thin, homogeneous, 
isotropic, of constant thickness over the whole tube length with linear distensibility 
and negligible compressibility to ensure that the same characteristic mechanical 
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response applies to all cross sections at all axial positions and in all radial directions 
leading to an axi-symmetric expansion in the tube geometry. We also rule out the 
possibility of a collapsed tube by imposing the condition that the pressure-modified 
cross sectional area cannot be smaller than its reference state corresponding to 
the reference pressure where the tube assumes its regular cylindrical shape under 
infinitesimal stress. 

Although the individual tubes in the network can in principle have different 
pressure-area mechanical characteristics, for the sake of simplicity the results gen- 
erated and reported in this paper are obtained with networks in which all the tubes 
are assumed to have the same mechanical characteristics. This is normally the case 
in the most common types of synthetic flow networks. It is also realistic in most bi- 
ological networks, such as blood circulation system, as long as the modeled network 
does not span different categories of vessels belonging to different circulation sub- 
systems with different mechanical characteristics, e.g. arteries and veins. Anyway, 
the model can be easily modified to accommodate multiple tube wall mechanical 
characteristics if such a modification is required. 

As done in the past by many network modelers [10-14], the network models can 
also be used as prototypes for porous media. This can be very good approximation 
when the void space is characterized by rather regular patterns with statistical 
distributions that can be well averaged and described by the distribution of the 
tubes in the network. 


2 Method 

According to the residual based pore scale network modeling approach, the flow rate 
and pressure fields in a fully connected network of tubes with certain mechanical 
properties that is filled with a particular fluid and subjected to a pressure gradient 
can be obtained by imposing the presumed boundary conditions on the inlet and 
outlet boundary nodes and enforcing a flow continuity condition on each interior 
node of the network using a characteristic flow relation. This relation normally 
correlates the volumetric flow rate to the boundary pressures in a single conduit and 
is based on the particular rheological fluid model and the characterizing properties 
of the flow conduit such as the flow of a power law fluid in a tube of elastic nature. 
Since in most cases these problems are non-linear requiring the solution of a system 
of nonlinear simultaneous equations, a non-linear iterative solution scheme, such 
as Newton-Raphson method, is then employed by starting from an initial guess for 
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the flow and pressure fields and iterating, while trying in each iteration to minimize 
the residual flow at each interior node by imposing the continuity condition. The 
final solution is then obtained when two consecutive iterations converge to the same 
solution within a given error tolerance. This method is fully described in [9] for 
the case of Newtonian flow in distensible networks. Further technical details about 
the non-linear residual based solution scheme can be found in [15], although this 
reference is related to a one-dimensional Navier-Stokes finite element flow model 
rather than a pore scale network model. 

Regarding the analytical relations correlating the pressure to the volumetric 
flow rate in a single tube that to be used in the proposed pore scale model for 
the flow of power law fluids in elastic networks, they are given by the following 
equations as derived previously in [8] 
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The first of these equations is based on a mechanical elastic response for the tube 
given by 


P = 7 (A - A 0 ) 

while the second is based on an elastic tube model given by 


( 3 ) 



( 4 ) 


In the last four equations, Q is the volumetric flow rate, 7 and /3 are the stiffness 
proportionality factors, A is the tube cross sectional area corresponding to the 
actual pressure p, A 0 is the reference cross sectional area corresponding to the 
reference pressure which for simplicity is assumed to be zero, p in and p ou are the 
pressure at the tube inlet and outlet respectively, L is the tube length, and k and 
n are the power law consistency factor and flow behavior index respectively. 
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3 Implementation and Validation 


The network flow model described in this paper was implemented in a pore scale 
network modeling computer code which is based on a Newton-Raphson iterative 
solution scheme that employs several numerical solvers (UMFPACK, LAPACK, 
SPARSE, and SUPERLU). An extensive series of tests were then carried out to 
inspect and assess the model. In the following we summarize some of the results 
of these tests in the context of model validation. As we have no experimental data 
to compare with, the validation tests are mainly focused on observing qualitatively 
and quantitatively correct trends related to special and limiting cases where the 
right outcome can be predicted in advance. As most of the validation tests rely on 
pore scale network flow models that are based on different underlying characteristic 
flow equations, we outline these equations in the following paragraph. 

In reference [8], the following two Equations for the flow of Newtonian fluids in 
elastic tubes, whose derivation is based on the same method and assumptions used 
in the derivation of the power law formulae (i.e. Equations 1 and 2), were derived 

0 _ (j Pin + lApf - ( p ou + yAp) 3 
247r/i7 2 L 

and 


Q 


P 

407r/iA o L 


A o \ 1 

^ Pin + \/A I 



(6) 


where p is the fluid dynamic viscosity. These equations are based on Equations 3 
and 4 respectively. In our validation tests we used these two formulae that describe 
the characteristic flow in the network tubes and whose derivation is completely 
independent from the derivation of the power law formulae to test the behavior 
of the power law flow model in this special case since the Newtonian fluids are a 
sub-category of the power law fluids corresponding to n = 1. We also used the 
following equation, which correlates the flow rate to the pressure drop for the flow 
of power law fluids in rigid cylindrical tubes [14, 16], to characterize the flow in the 
network tubes in some of these tests 


I / , \ 3+1 /n 
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where x is the tube axial coordinate. Finally, the well known Poiscuille model was 
also used in one of our validation tests. 


5 


The following points outline the correct trends that have been observed and 
hence been used to validate the model where all these trends are verified for both 
p-A elastic tube models as given by Equations 3 and 4 

1. The convergence of the solutions of the power law elastic network model 
(characterized by Equations 1 and 2) to the solutions of the power law rigid 
network model (characterized by Equation 7) by increasing the tube wall 
stiffness. These trends are demonstrated in Figure 1 using a computer gen- 
erated fractal network. Full description of this type of networks is provided 
in [9]. 

2. The convergence of the solutions of the power law elastic network model 
(characterized by Equations 1 and 2) to the solutions of the Newtonian clas- 
tic model (characterized by Equations 5 and 6 respectively) when n — 1. 
These trends are demonstrated in Figure 2 using a computer generated or- 
thorhombic network. Full description of this type of networks is provided in 
[9], 

3. As a consequence of the previous two points, it is expected that the solutions 
of the power law elastic network model would converge to the solutions of the 
Newtonian rigid network model (Poiscuille) when n = 1 and the stiffness of 
the tubes wall is high. This trend is also verified as demonstrated in Figure 
3 using a computer generated fractal network. 

4. Another validation test is the agreement between the solutions of a discretized 
tube, as obtained numerically from the residual based pore scale network 
modeling approach, and the analytical solutions for a single tube, as given 
by Equations 1 and 2. This is verified and demonstrated in Figure 4. A 
similar test corresponding to the Newtonian case has confirmed the agreement 
between the numerical solutions of the Newtonian elastic network model and 
the single tube analytical expressions as given by Equations 5 and 6. 

4 Assessment and Comparison 

It is difficult to make a fair comparison between two mathematical or computational 
models that are based on different theoretical and numerical frameworks. This 
applies to any comparison between our proposed pore scale network model and any 
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Figure 1: Convergence of solutions of power law pore scale elastic network model 
with high tube wall stiffness (using Equations 1 and 2 for their single tube char- 
acteristic flow) to their equivalents of rigid model (using Equation 7 for single 
tube characteristic flow) (a) for the first p-A elastic model (Equation 3) with 
7 = 10 12 Pa.m -2 using a power law fluid with n = 0.7 and k = 0.025 Pa.s n ; 
and (b) for the second p-A elastic model (Equation 4) with /3 = 10 6 Pa.m using a 
power law fluid with n = 0.8 and k = 0.015 Pa.s n . In both cases a fractal network 
consisting of 511 tubes in 9 generations with an area-preserving branching expo- 
nent of 2 [9, 17] and an inlet main branch with L = 0.08 m and R 0 = 0.016 m was 
used. The vertical axis in the sub-figures represents the net inflow/outflow, Q, in 
m 3 .s -1 while the horizontal axis represents the inlet boundary pressure, p in , in Pa. 
The outlet pressure in both cases is zero. 


one of the existing models such as the one-dimensional Navier-Stokes flow model 
which is usually associated with a finite element numerical implementation [15] . In 
reference [8], we outlined some issues related to such a comparison with regard to 
the flow in single tubes. Most of that discussion also applies to the network flow 
models that are based on these two distinctive formulations. 

However, there are some obvious advantages in using the pore scale network 
modeling approach over the one-dimensional Navier-Stokes finite element model. 
Some of these advantages are based on our personal experience and hence may be 
related to our computational frameworks for these models and their particular im- 
plementations. The outcome may be different for other implementations although 
we believe that the assessment is valid in general. 

A major advantage of the proposed pore scale network modeling approach is 
the incorporation of non-Newtonian effects, whereas the one-dimensional Navier- 
Stokes model is limited to the Newtonian flow. Also, since the pore scale modeling 
approach is based on analytically derived characteristic flow relations, it produces 
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Figure 2: Convergence of solutions of power law elastic network model (using 
Equations 1 and 2 for their single tube characteristic flow) to their equivalents 
of Newtonian elastic model (using Equations 5 and 6 for their single tube char- 
acteristic flow) when n — 1 (a) for the first p-A elastic model (Equation 3) with 
7 = 10' Pa.m -2 using a fluid model with k — n — 0-03 Pa.s; and (b) for the 
second p-A elastic model (Equation 4) with /3 = 80 Pa.m using a fluid model with 
k = p = 0.008 Pa.s. In both cases an orthorhombic network consisting of 946 tubes 
with an average R a of 0.0042 m and a standard deviation of 0.0015 m [9] was used. 
The axes and outlet pressure are as in Figure 1. 


exact solutions within its defining assumptions considering that the allowed error 
tolerance in the residual based iterative scheme can be minimized to the limit of 
machine precision and hence can be forced to vanish for all practical purposes. 
Computationally, we observe that the proposed pore scale modeling approach is 
superior in terms of robustness, reliability, ease of implementation, low compu- 
tational overhead, and better rate and speed of convergence. For example, we 
experienced instabilities and convergence difficulties [18] when using the first p-A 
relation with the one- dimensional Navier-Stokes finite element flow model but not 
with the power law pore scale network model even though the opposite should be 
expected since the latter employs a non-linear non-Newtonian rheology while the 
former employs a linear Newtonian rheology. Similar assessment should also ap- 
ply in general to comparisons between the proposed pore scale network model and 
other existing models. 


5 Conclusions 

The flow of non-Newtonian fluids in distensible networks and porous media is 
investigated where a residual based pore scale network modeling method employing 
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(b) Second p-A model 


Figure 3: Convergence of solutions of power law pore scale elastic network model 
with high tube wall stiffness and n = 1 (using Equations 1 and 2 for their single 
tube characteristic flow) to their equivalents of rigid Poiseuille model (a) for the 
first p-A elastic model (Equation 3) with 7 = 10 11 Pa.m' 2 using a fluid model with 
k — fi — 0.005 Pa.s; and (b) for the second p-A elastic model (Equation 4) with 
/3 = 10 6 Pa.m using a fluid model with k = p = 0.003 Pa.s. In both cases a fractal 
network consisting of 255 tubes in 8 generations with a branching exponent of 2.5 
[9, 17] and an inlet main branch with L = 0.07 m and R 0 = 0.0105 m was used. 
The axes and outlet pressure are as in Figure 1. 



Figure 4: Comparison between the analytical solutions, as given by Equations 1 
and 2, for a single elastic tube with L = 0.6 m and R 0 = 0.08 m, and the numerical 
network model solutions for the same tube obtained by discretizing the tube in the 
axial direction into 101 rings and treating the discretized ensemble as a network 
of serially connected tubes (a) using the first p-A elastic model (Equation 3) with 
7 = 50000 Pa.m -2 , n = 0.8, and k = 0.2 Pa.s n ; and (b) using the second p-A 
elastic model (Equation 4) with (3 = 7500 Pa.m, n = 1.2, and k = 0.08 Pa.s n . The 
axes and outlet pressure are as in Figure 1. 
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a recently derived formulae for the flow of power law fluids in clastic cylindrical 
tubes is proposed. This modeling method was implemented in a computer code and 
initial results based on extensive tests were obtained. The results show a number 
of sensible trends which include convergence to special and limiting cases where 
the outcome can be predicted from previously validated alternative models. 

The investigation also included assessment of the proposed model and general 
comparison to the existing ones notably the finite element based one-dimensional 
Navier-Stokes model. Several advantages of the proposed model have been asserted 
such as producing exact solutions within its underlying assumptions and allowed 
numerical errors, reliability, ease of implementation, low computational cost, and 
fast and smooth convergence. 

Although there are previous attempts where certain non-Newtonian effects were 
considered in the flow through small scale distensible branching flow ensembles, the 
main novelty of the proposed method is that it considers non-Newtonian rheology 
directly and on a large network scale employing analytically derived underlying 
flow equations. 

Nomenclature 

f3 stiffness factor in one of the pressure-area elastic models 
7 stiffness factor in another pressure-area clastic model 
p fluid dynamic viscosity 

A tube cross sectional area 

A a tube reference cross sectional area at reference pressure 
k consistency factor of power law fluid 
L length of tube 

n flow behavior index of power law fluid 
p axial pressure 

Pin inlet pressure 

p ou outlet pressure 
Q volumetric flow rate 

R 0 radius of tube corresponding to tube reference area 
x tube axial coordinate 
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